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Abstract 


A procedure has been developed for modeling wind tunnel flows using computa- 
tional fluid dynamics. Using this method, a numerical study has been undertaken 
to explore the effects of solid wind tunnel wall proximity and Reynolds number on a 
two-dimensional airfoil model at low speed. Wind tunnel walls are located at varying 
wind tunnel height to airfoil chord ratios and the results are compared with freest ream 
flow m the absence of wind tunnel walls. Discrepancies between the constrained and 
unconstrained flows can be attributed to the presence of the walls. Results are for 
a Mach Number of 0.25 at angles of attack through stall. A typical wind tunnel 
Reynolds number of 1,200,000 and full-scale flight Reynolds number of 6,000,000 
were investigated. At this low Mach number, wind tunnel wall corrections to Mach 
number and angle of attack are supported. Reynolds number effects are seen to be a 
consideration in wind tunnel testing and wall interference correction methods. 

The study uses an unstructured grid Navier-Stokes code with Baldwin- Lomax 
turbulence model. The numerical method is described since unstructured flow solvers 
present several difficulties and fundamental differences from structured grid codes, 
especially in the area of turbulence modeling and grid generation. 


ii 


PRECEDE FAGZ CLANK NOT FILMED 




Acknowledgements 


This work was performed under Boeing contract Y429804-0957N and NASA contract 
NCC2-055. Their support is gratefully acknowledged. I would like to thank Tim 
Barth and Marshad Merriam at NAS A- Ames for the use of their flow solver and 
grid generation codes as well as numerous helpful discussions. 




Contents 


Abstract iii 

Acknowledgements iv 

Contents v 

List of Tables vii 

List of Figures viii 

Nomenclature x 

1 Introduction 1 

1.1 Applications 1 

1.2 Current Configuration 2 

2 Computational Method 5 

2.1 Equations 5 

2.2 Discretization 6 

2.3 Fluxes 7 

2.4 Higher Order Accuracy 9 

2.5 Turbulence Modeling 11 

2.5.1 Baldwin- Lomax Model 11 

2.5.2 Implementation on Unstructured Grids 12 

2.5.3 Considerations for Complex Flows 13 


iv 




2.6 Time Evolution 


14 


3 Grid Generation 10 

3.1 Method 16 

3.2 Grid Adaption 25 

3.3 Grids 27 

4 Boundary Conditions 32 

4.1 Solid Walls 32 

4.2 Inflow and Outflow 34 

5 Discussion 36 

5.1 Results 36 

5.2 Lift 4g 

5.3 Drag 54 

5.4 Code Characteristics 56 

6 Conclusions 58 

References 61 


v 




List of Tables 

1 Grid Dimensions: h/c — 2.25, Re = 1,200,000 27 

2 Lift and Drag Coefficients 51 

3 Pressure and Viscous Drag Coefficients 54 


vi 




List of Figures 

1 Mesh Geometry and Control Volumes 7 

2 Control Volume and Variables g 

3 Higher Order Accuracy, Nodal Gradient Calculation 10 

4 Overlaid Structured Grids yj 

5 Overlapping Sections Removed from Structured Grids L8 

6 Filtered Points and Spacing Function 20 

7 Filtered Points in the Wake 21 

8 Delauney Triangulation 22 

9 Non-Optimal Delauney Triangulation with Viscous Grids 23 

10 Edge Swapping in Locally Stretched Plane 24 

11 Final Grid - Quadrilaterals and Triangles 26 

12 Grid Adaption 28 

13 Grids for h/c = 2.25, Re = 1,200,000 29 

14 Grid for Unconstrained Flow 31 

15 Weak Enforcement of Boundary Conditions 33 

16 Airfoil Pressures: a = 0 degrees, Re = 1,200,000 37 

17 Airfoil Pressures: a = 0 degrees 37 

18 Airfoil Skin Friction: a = 0 degrees, Re = 1,200,000 38 

19 Airfoil Skin Friction: a = 0 degrees, Re = 6,000,000 38 

20 Wall Pressures: a = 0 degrees 39 

21 Airfoil Pressures: a = 5 degrees, Re = 1,200,000 40 

22 Airfoil Pressures: a = 5 degrees 41 

23 Airfoil Skin Friction: a = 5 degrees, Re = 1,200,000 42 


Vll 




24 Airfoil Skin Friction: a = 5 degrees, Re = 6,000,000 42 

25 Wall Pressures: a = 5 degrees 43 

26 Mach Number Contours: a = 5 degrees, hjc = 1.5, Re = 1,200,000 . . 44 

27 Mach Number Contours: a = 5 degrees, unconstrained, Re = 1,200,000 44 

28 Airfoil Pressures: a = 10 degrees, Re = 1,200,000 45 

29 Airfoil Pressures: a = 10 degrees 45 

30 Airfoil Skin Friction: a = 10 degrees, Re = 1,200,000 46 

31 Airfoil Skin Friction: a = 10 degrees, Re = 6,000,000 46 

32 Wall Pressures: a = 10 degrees 47 

33 Airfoil Pressures: a = 15 degrees, Re = 1,200,000 47 

34 Airfoil Pressures: a = 15 degrees 48 

35 Blow Up of Leading Edge Mach Number Contours: a = 15 degrees, 

hjc = 1.5, Re = 6,000,000 49 

36 Airfoil Skin Friction: a = 15 degrees, Re = 1,200,000 50 

37 Airfoil Skin Friction: a = 15 degrees, Re = 6,000,000 50 

38 Wall Pressures: a = 15 degrees 51 

39 Mach Number Contours: a = 15 degrees, hjc = 1.5, Re — 1,200,000 . 52 

40 Mach Number Contours: a = 15 degrees, unconstrained, Re = 1,200,000 52 

41 Lift Curves versus hjc 53 


Vlll 




Nomenclature 


A 

c 


C d 

Cdr,.. 

C^. c 

Cf 

Ct 

c p 


F, G 

hjc 


j 

k 

l 

M 


n 

n,t 

P 

PE 

Q 

P-u P2 


Area 

airfoil chord, 8 in. 

speed of sound 

drag coefficient, j — ^-5— 

2 P°° U oo C 

pressure or form drag coefficient 
viscous drag coefficient 
skin friction coefficient, » T » “ a “ 
lift coefficient, ■» — ^-=— 
pressure coefficient, t — *-?- 
energy 

Cartesian flux vectors 

total wind tunnel height to airfoil chord ratio 
cell index 
node index 

coefficient of thermal conductivity 
turbulent mixing length 
Mach number 
time index 

normal, tangential coordinates 
static pressure 
potential energy 

vector of conserved variables (mass, momentum, energy) 
Riemann invariants 


xi 




R residual 

Re Reynolds number, p^u^c/ 

s entropy 

5 length 

t time 

T static temperature 

u y v Cartesian velocity components 

friction velocity, yjT wall / p wall 
x,y Cartesian coordinates 

y + turbulent law of the wall nondimensional normal coordinate, u r y/i/ wall 


a 

7 

0 

\ 

P 

Pt 

v 

P 

T 

V 

Q 


angle of attack, degrees 
ratio of specific heats, 1.4 
angle between adjacent edges 
second coefficient of viscosity, —2/3/i 
coefficient of viscosity 
turbulent eddy viscosity 
kinematic viscosity, pj p 
density 

viscous shear stress 

vorticity, g - g 

control volume with boundary SQ 


Subscripts 

e Euler, inviscid 

I<,R states on either side of an edge (left, right) 

n,t normal, tangential components 

v viscous 

x,y Cartesian derivatives, components 

oo freestream value 


x 




Chapter 1 
Introduction 


This work describes a procedure which has been developed for modeling wind tunnel 
flows using computational fluid dynamics (CFD). Possible applications of the CFD 
code, grid generation programs, and boundary conditions discussed here are numer- 
ous. Results of a fundamental study of the effects of solid wind tunnel wall proximity 
and Reynolds number on an airfoil model are presented. 


1.1 Applications 

CFD validation benefits from the ability to numerically calculate wind tunnel flows. 
When wind tunnel walls are simulated directly, the errors due to wind tunnel wall 
interference correction methods are eliminated. This enables a direct comparison of 
experimental data with computational results. With the effects of the wall correction 
method removed, discrepancies can more accurately be attributed to numerical errors. 
This idea is promising for the validation of turbulence models in that it can pinpoint 
areas of weakness in the model. Errors due to the simulation of porous or slotted 
walls can be eliminated if solid walls are used. 

A broad range of flow regimes can be handled more expediently with a CFD code 
than with experimental tests. The numerical simulation also provides considerably 
more data than is available from wind tunnel tests. Alternatively, CFD simulations 
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can be used for the development of wind tunnel wall interference correction meth- 
ods. The computational data can determine the effects which must be included in a 
correction method and the appropriateness of a particular method. Although wind 
tunnel wall interference assessment /correction (WIAC) techniques have made consid- 
erable progress and are often quite reliable, numerous issues must still be addressed 
in future developments [1, 2]. Reynolds number effects are an important area which 
can be investigated using a Navier-Stokes code. For CFD validation methods and 
wind tunnel wall interference correction methods which employ wind tunnel mea- 
surements directly, the large amount of data available from CFD analyses can be 
used to determine the amount and location of experimental data that is required so 
that experimental measurements will be minimized. CFD codes are often used as 
part of WIAC methods to determine effective airfoil shapes and to iterate on Mach 
number and angle of attack corrections [1]. 

Although the eventual aim of a wind tunnel experiment is to be able to relate it 
to full-scale free flight, understanding of both flow conditions is necessary. The CFD 
code and grid generation system developed here allow for studies of Reynolds number 
effects; solid, porous, and slotted wall effects; wall proximity effects; laminar versus 
turbulent effects; as well as the modeling of a wide range of geometries, angles of 
attack, and flow conditions. 


1.2 Current Configuration 

The code was used in this work to study the effects of wind tunnel wall proximity 
and Reynolds number on a two-dimensional single element airfoil model at low speed. 
Wind tunnel wall proximity effects are important because it is desirable to put as 
large a model in a wind tunnel as possible. Conversely, small wind tunnels are more 
economical to operate. Larger models are not only easier to build and instrument, 
but they more closely match full-scale Reynolds numbers. Reynolds number and wall 
proximity effects are coupled in that attempts to duplicate higher Reynolds numbers 
are often limited by blockage effects due to wall proximity. 

The computations in this work use a Boeing advanced transport research airfoil. 
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The airfoil has camber. Although the model includes a main airfoil and flap, the flap 
has been retracted for this study. Four angles of attack were computed: a = 0, 5, 10, 
and 15 degrees. The first three cases yield attached flow while at fifteen degrees the 
airfoil is mostly separated. 

Solid wall wind tunnel flows were computed with total height to airfoil chord 
ratios of 1.5, 2.25, 4.5 and oo (unconstrained). The unconstrained case sets the outer 
boundary 12 chords away from the airfoil. For conformity with the wind tunnel cases, 
the grid for the unconstrained case can be equated to a height to chord ratio of 24, 
although the boundary conditions are different. The wind tunnel inflow plane was 
located five chords upstream from the leading edge and the exit is fifteen chords 
downstream from the leading edge. The wind tunnel walls were modeled as solid. 

The flow is at a freestream Mach number of 0.25. Two Reynolds numbers based 
on the airfoil chord were considered: 1,200,000 and 6,000,000. A Reynolds number of 
1,200,000 is typical of a low speed wind tunnel as is a height to chord ratio of 2.25. 
This configuration is patterned after experimental work done on the airfoil with flap 
extended at the Stanford subsonic wind tunnel [3]. Full-scale testing is simulated by 
the higher Reynolds number. 

The CFD code is an unstructured grid Navier-Stokes solver with a Baldwin- 
Lomax algebraic turbulence model. The Reynolds-averaged Navier-Stokes equations 
are solved so that viscous phenomenon, most importantly the wake and boundary 
layer, can be included and investigated. Reynolds number effects have been shown to 
affect airfoil characteristics and the level of wall interference [2]. A turbulence model 
is used because the flows studied here are typically turbulent. In the wind tunnel, 
transition occurs naturally or the boundary layer is artificially tripped to simulate 
full-scale conditions. 

A major advantage of an unstructured grid is its flexibility in modeling complex 
geometries. The single element airfoil investigated here is a precursor to future work 
on a multielement model for which an unstructured grid is well suited, although tur- 
bulence modeling can become more complex. The grid generation scheme is demon- 
strated using a multielement airfoil and wind tunnel configuration. Another benefit 
of using unstructured grids is the ease with which they can be refined and adapted 
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to resolve flow features. Currently, to obtain the benefits of using unstructured grids, 
while removing the penalties of structured grids, extra steps are required in the grid 
generation process. The computational method will be discussed along with the ap- 
propriate boundary conditions for wind tunnel flows, particularly for the inflow and 
outflow planes. 



Chapter 2 

Computational Method 


The computational fluid dynamics code solves the Navier-Stokes equations in con- 
servative form on an unstructured grid in two dimensions. The space derivatives are 
discretized using a finite volume formulation which is employed up to the boundaries. 
Roe’s flux difference splitting is used for the inviscid terms and central differencing is 
used for the viscous terms. The Baldwin-Lomax algebraic turbulence model is imple- 
mented. Because of the unstructured grid, explicit time stepping is used to advance 
the solution to steady state. Several enhancements are added to speed convergence. 
The code was developed by Timothy Barth at NASA-Ames Research Center. See 
Reference 4 for details and further references. 


2.1 Equations 

The Navier-Stokes equations in differential conservative form are 

dQ «3F e <9G e _ 5Fv <9G V 

dt + dx dy dx dy 

where, 

Q = [p,pu,pv,e] T 
F e = \pu,pu 2 +p,/mt;,u(e + p)] 

G e = {pv^puv^pv 2 + p,v(e + p)] 


( 1 ) 
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Txy Tyx 


Fy — [0 j Txx * Try y ^^"xx "t” ^^*xy ~t“ ^-^x] 

G v = [0, ^yx ? ^yy ? ^^"yx H” ^^yy 
T xx = 2 \LU X - A(u x + v y ) 

j. l^Uy "f* Vx ) 

Tyy - 2/iVy ~ A(ti* + Vy) 


In integral form they are written 

dQ 


L 


<9F e dG , 

+ 


L ) 


dFy OGy 


+ 


dt dA + Jn\'fc"~d71""~Jn\dz ' 8y 


dA 


( 2 ) 


where Q represents a control volume with differential area dA. The area integrals can 
be converted to surface integrals over the boundary of f2, resulting in a flux balance 
over the control volume: 

— f Q dA + [ (n x F e + n v G e ) dS = f (n x F v + n y G v ) dS (3) 

dt Jo Jao Jan 

where dS is a differential length on the boundary dQ, and n x and n y are the unit 
normals to the boundary. 


2.2 Discretization 

The finite volume formulation used here discretizes the control volume as a polygon 
so that the integration becomes a numerical quadrature over the discrete faces of the 
polygon: 


d Qj 

dt 


nfacet(j) 


+ Aj x J2 MF«-Fv) + ny(G.-G T )]A5 i = 0 


(4) 


t— 1 


Arbitrary polygons are handled in this manner. The term in brackets in Equation 4 
is the oriented flux along the i th edge of the control volume with length AS,. The 
vector of conserved variables Q at node j is the cell averaged data such that 


= -v L q 


dA 


(5) 


The flow variables are, therefore, piecewise constant in each cell. Data is stored at 
the vertices of the mesh, and the median dual mesh is used as the control volumes. 
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mesh 

centroid dual 
median dual 



Figure 1: Mesh Geometry and Control Volumes 

The median dual is formed by connecting cell centroids to edge centers as shown in 
Figure 1. Currently, the grid cells are either triangles or quadrilaterals. The numerical 
integration is done over all discrete segments of the control volume for better accuracy 
at triangle/quadrilateral cell interfaces. This is necessary because the difference in 
truncation error between the two cell geometries causes accuracy problems. This 
improved quadrature adds additional computational time. 

2.3 Fluxes 

The inviscid fluxes are computed using Roe’s flux difference splitting. Rewriting 
Equation 4, the oriented flux, f n , is considered to be a function of the states on either 
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control volume 



Figure 2: Control Volume and Variables 

side of an edge of the control volume, Ql and Qr (Figure 2): 

JA, nfacf(j) 

+ Aj 1 £ f n (Qt,QR)AS, = 0 (6) 

at i =1 

Roe’s flux function defines the oriented flux to be an approximate solution to the 
Riemann problem at this edge: 

fn(Qlo Qr) = 2 [^n(Ql«) + ^(Qtt)l — ^ Q®-)l (Q L - ^ 

The first term represents a standard central difference while the second term adds an 
upwind influence by distinguishing incoming and outgoing waves. This term also has 
the effect of adding artificial dissipation to stabilize the solution. The matrix A is the 
Jacobian matrix of f, dl/dQ. The matrix |j 4| has the same eigenvectors as A, but its 
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eigenvalues are the absolute value of those of A. Qj, and Qr determine the order of 
accuracy of the scheme and are defined in the next section. 

In the viscous fluxes the values of the conservative variables at an edge of the 
control volume are an average of the nodal values on either side of the edge. The 
gradients in the viscous fluxes there are computed similar to a finite difference formu- 
lation using the values Qji and Qj 2 (Figure 2). In addition to the gradient component 
along the jl~j2 edge, the perpendicular component of the gradient in cell i is added: 

V Q v = Qj *~ t Qjl n + VQi.t (8) 

2.4 Higher Order Accuracy 

A first order scheme is created if Ql and Q& in the inviscid fluxes are taken as 
nodal values. A higher order scheme is created if the cell averaged nodal data is 
reconstructed to be piecewise linear in a cell rather than assumed constant. In the 
piecewise linear .case, Ql and Qr at an edge of the control volume are determined 
from an expansion about nodal data: 

Ql = Qj + VQj • Ar (9) 

where Qj and VQj are the nodal data and gradients and Ar is the vector between 
the node and the midpoint of an edge of the control volume as seen in Figure 3. The 
nodal gradients used in reconstruction are computed from: 

/ VQ dA = <£ Qhds (10) 

Jo Jan 

With linear reconstruction VQ is constant so that 

VQ; = Aj l <£ Qn ds (11) 

J Jdu 

Using a path of integration made up of the surrounding cells as in Figure 3, the nodal 
gradient is computed as 

nfacet(j) 

VQj = ^7 l £ Qi nj As (12) 

i= 1 

The cell gradients, VQj, used in the viscous fluxes are computed in a similar manner 
using the cell boundary as the contour of integration. 
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path of integration for 
nodal gradients 



Figure 3: Higher Order Accuracy, Nodal Gradient Calculation 
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2.5 Turbulence Modeling 

Barth has implemented the Baldwin-Lomax algebraic turbulence model [5] in the 
code. The effects of turbulence are modeled by the addition of a turbulent eddy 
viscosity to the molecular viscosity in the viscous stresses. This requires special 
procedures to be executed on an unstructured grid. 

2.5.1 Baldwin-Lomax Model 

The two layer model is defined as follows. In the inner layer: 

Pt inn „ = />* 2 M ( 13 ) 


where, 

( = ky[ l-e~ y+/A+ ] 

is the turbulent mixing length. In the outer layer: 

. = KC C PpF wakt F Kl ' b {y) ( 14 ) 


where, 

F\uakt = min(?/max F max , CwKVmaxU j F max ) 
F{y) = y\w\[l- e~ v+/A *} 



Ufa q — w maje n m j n 


The constants k, >1 + , K, Cep , Cwk > and Ckia are those defined in the original 
paper by Baldwin and Lomax [5]. The coordinate direction y is actually the normal 
distance from the wall and is not always aligned with the Cartesian direction. The 
crossover point from inner to outer layer is determined by the location where Pt, nn . r 
exceeds in a profile. The quantity F is the damped moment of vorticity, and 

the location of its maximum in a boundary layer or wake profile is used as a length 
scale to compute turbulent viscosity. In structured meshes it is a simple procedure to 
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scan boundary layer or wake profiles to find this maximum because the profiles are 
typically coordinate lines normal to the body or wake centerline. In an unstructured 
grid, these lines may exist, especially if the unstructured grid was generated from a 
structured one, but they are not easily accessed since unstructured meshes have no 
preferred set of coordinate directions. 

2.5.2 Implementation on Unstructured Grids 

The current implementation constructs these lines using geometry information from 
only the boundary edges and points. The procedure for near wall flows is as follows: 
For each grid point find the nearest no slip boundary edge, the distance to the nearest 
point on this boundary edge, and an interpolation factor based on the location of this 
nearest point on the boundary edge. The interpolation factor is used to determine 
which boundary edge endpoint is referenced and to interpolate to the normal line 
represented by this endpoint. The trailing edge point is marked and all points which 
have this point as the closest are flagged as wake points. A mixing length model is 
used in the wake. All points which reference a particular boundary point as closest 
are considered to be in that point’s boundary layer profile. The points are sorted in 
order of increasing distance from the body in a preprocessing stage. 

At each time step the moment of vorticity, F, and the inner layer turbulent vis- 
cosity are computed for all points. The points are then processed in sorted order to 
find the maximum moment of vorticity, F ma xi f° r each boundary layer profile and, 
therefore, each boundary point. This is possible since the interior points where the 
function is computed all reference a boundary point. The value of the function in 
the field is used in conjunction with the interpolation factor to compute properties at 
specific boundary stations. The outer layer turbulent viscosity for each interior point 
is calculated as a function of the distance of the point from the wall and F max of the 
referenced boundary point. Finally, the points are again processed in sorted order 
to find the crossover point from inner to outer layer for each profile. No transition 
model is implemented, but this is not a limitation of the unstructured grid. 

The method works particularly well when the unstructured mesh is originally 
obtained from a structured mesh and lines normal to the body do exist. If the grid is 
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completely unstructured the current implementation can break down. This happens 
when a boundary point cannot generate a complete set of points which span the 
entire profile. An incorrect length scale can be computed because the boundary point 
does not see the outer edge of the layer. A more robust implementation would have 
boundary points pointing to field points to make up a profile rather than field points 
pointing to boundary points. Rostand [6] does this by generating normal lines at each 
boundary point and storing the interior edges and points which the profile intersects. 
All interior points then use boundary normal lines to interpolate from. An alternative 
method for creating boundary normals is to employ background structured grids and 
interpolate flow variables and turbulence quantities back and forth, but this can be 
quite memory intensive [7]. No problems were encountered with the current method 
since the grids had some original structure to them, and when boundary points were 
added, care was taken to add a complete profile to go with them. 

The process of setting up the nearest boundary edges and points as well as the 
sorting is only performed once because this information is stored. Storage required 
for the turbulence model is two real arrays: 1) the distance to the nearest wall and 
2) the interpolation factor, and two integer arrays: 1) the referenced boundary edge 
and 2) the order of sorted points. All of these arrays have dimensions of the number 
of nodes in the mesh. These arrays are the only additional storage requirements 
due to the unstructured nature of the grid. The Baldwin-Lomax model requires the 
storage of the maximum moment of vorticity and its location as well as inner and 
outer turbulent viscosities, but these arrays represent no additional storage because 
previously dimensioned arrays are reused for these variables. The current method 
uses less than 5% of the total memory requirements. 

2.5.3 Considerations for Complex Flows 

Confluent or overlying shear layers which are present in multielement airfoil flowfields 
are currently not handled correctly. When using an algebraic model, the boundary 
layer of the flap must have a different length scale than the overlying main element 
wake. This is not accounted for as the the current implementation only creates one 
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length scale per profile. If an algebraic model is used several alternatives can be pro- 
posed to treat the interaction of turbulent wakes and boundary layers. Mavriplis [7] 
has implemented the aforementioned method of background structured grids and has 
shown impressive results for a four element airfoil. Alternatively, with minor modi- 
fications to the current method, wake edges can be included in the list of boundary 
edges. The turbulence model will then perform in a manner similar to structured grid 
implementations which use a grid line approximating the wake centerline. When adja- 
cent boundary layers and wakes start to merge, a combination of turbulent viscosities 
is required. 

For complex flows such as confluent wakes, algebraic turbulence models are prob- 
ably no longer reliable. Execution of a one or two-equation model such as k-e presents 
little difficulty on an unstructured mesh. One of the major disadvantages is the time 
required to integrate another, possibly stiff, equation. A modified k-t model has 
performed reasonably for merging and confluent wakes and boundary layers [8]. 


2.6 Time Evolution 


The time derivative in Equation 6 is discretized using an explicit 3 stage Runge-Kutta 
time stepping scheme. For the finite volume space discretization: 


Qj (0) = 

Qj (n) 

dt maXj 
Cl V olj 


Qj (1) = 

Qj (0) - 


Qj (2) = 

Qj (0) - 

~ dtfnaxj 

° 2 V olj 

■Rj(Qj (1) ) 

Qj (3) = 

Qj (0) - 

/-* dt maXi 
° 3 Volj 


Qj (n+1) = 

Qj (3) 



Ci = .18, 

C 2 -- 

= .5, 

C 3 = 1.0 


To accelerate convergence local time stepping is used. The time step in each cell 
is determined by the local CFL number in that cell and, therefore, varies throughout 
the mesh. The calculation is no longer time accurate. The CFL number is computed 
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from monotonicity principles and is typically overrestrictive when compared with a 
conventional CFL number [4]. 

Another method used to speed convergence is grid sequencing. The solution is 
initially started on a coarse grid where the Euler equations are solved. This allows the 
solution which uses freestream initial conditions to set up rapidly. The viscous terms 
are then turned on to set up a boundary layer and wake. Grid adaption is used to 
refine the wake and boundary layer until the desired spacing is obtained. Examples 
of this procedure are given in the following section on grid generation 
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Grid Generation 


The major advantage in using unstructured grids is the inherent flexibility in treat- 
ing complicated geometries. Complex flows can also be computed more easily and 
accurately using grid adaption. Extra effort, though, is required to obtain the un- 
structured grids suitable for viscous calculations. 

3.1 Method 

The method will be demonstrated to generate a fine mesh about a multielement 
airfoil in a wind tunnel. This configuration would be a difficult task for a structured 
grid generator. The first step in generating the unstructured grid involves generating 
structured grids about all components: main airfoil, flap, and wind tunnel [9]. The 
background wind tunnel grid is Cartesian while grids about the airfoils are body 
conforming C-meshes, best at resolving boundary layers and wakes. The grids are 
overlaid (Figure 4) and overlapping sections, for example, those which fall within a 
body, are removed (Figure 5). At this point quadrilateral cells in boundary layers and 
wakes can also be removed for later reinsertion into the unstructured grid. These large 
aspect ratio cells can cause problems in latter steps in the grid generation process. 
This sequence of steps is done interactively on a graphics workstation. 

The connectivity of the points is thrown out, and only the coordinates of the points 
are saved. These points are then filtered. Points that were required in the initial grid 
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Figure 4: Overlaid Structured Grids 
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Figure 5: Overlapping Sections Removed from Structured Grids 
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due to its structured nature but which are not necessary for accuracy are discarded. 
For example, clustering at the trailing edge in a structured grid extends out into the 
field where it is not needed. By setting up a function in space which takes on the value 
of the minimum specified spacing between points, points are automatically filtered. 
Figure 6 illustrates the resulting cloud of points after such a filtering. The size of the 
crosses at the points indicates the minimum spacing between points. Note the points 
that have been removed around the airfoil when the boundary layer cells were saved. 
Points are also removed from the far downstream wake where they are not needed as 
the wake dies out (Figure 7). 

The next step involves connecting the points using a Delauney triangulation al- 
gorithm [10] (Figure 8). A Delauney triangulation is constructed by making three 
points into a triangle if their circumcircle contains no other points. Another prop- 
erty of Delauney triangulation is that it maximizes the minimum angle in the mesh. 
The resulting grid is isotropic with no directional bias. This is not always optimal, 
especially for grids with high aspect ratio cells as is typically required for viscous cal- 
culations. Problems such as seen in Figures 8 and 9 can result when the triangulation 
does not recover the initial structured mesh. This is one of the reasons quadrilateral 
boundary layer cells are removed early in the process. Another reason is to help 
prevent edges from “breaking through” across a body. 

In order to regain a directionally biased mesh, after the global triangulation is com- 
pleted, the mesh is retriangulated locally in a stretched plane using an edge swapping 
algorithm [11]. The original structured grids give the information about the trans- 
formed plane, namely a stretching vector with magnitude and direction. The vector 
is associated with the aspect ratio and major axis of surrounding cells at a point. The 
maxmin angle property of Delauney triangulation is utilized locally by considering 
two triangles with a shared edge. The orientation of the shared edge is determined 
by satisfying the maxmin angle criterion in the stretched plane. Several iterations 
through the mesh result in edges being swapped into the desired configuration (Fig- 
ure 10). 

The mesh is now made up only of triangles. The next step involves removing 
diagonals to form quadrilaterals out of two adjacent triangles [10]. This recovers most 
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Figure 6: Filtered Points and Spacing Function 
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Figure 7: Filtered Points in the Wake 
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Figure 8: Delauney Triangulation 
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Figure 9: Non-Optimal Delauney Triangulation with Viscous Grids 
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of the structured grid that was not destroyed in the filtering process. At this point 
the boundary layer quadrilaterals are also merged back in (Figure 11). Quadrilaterals 
are used for several reasons. They work well in shear layers where triangles with small 
height to base ratios, which have edges that are not aligned with the flow gradients, 
need to be avoided. They are also more accurate than triangles. The flow code 
is edge based in that most operations are done in loops over edges. Quadrilateral 
cells reduce the number of edges thereby increasing speed and decreasing storage. 
Triangles are generally only used when it is necessary to transition between different 
mesh geometries, airfoil elements, or areas of grid refinement. 

The last step in the mesh generation process, which can be seen in Figure 11, 
was to smooth the points. Near boundaries with high curvature, crossed grid lines 
can result if a simple Laplacian smoothing, i.e. averaging of surrounding coordinate 
values, is used. The procedure employed places torsion springs between adjacent 
edges of a vertex and approximately solves a local minimization of potential energy 
problem for each node with its surrounding cells. To prevent very small cells and 
large area variations, a term proportional to the inverse area of a cell is also added 
to the potential energy, which is computed as 


PE(x h yj) 


nvertex(j) 

E 




e 2 { Xj , yj ) 


K 

M x i> Vi). 


(16) 


where A0 is the angle between two adjacent edges, Ai is the enclosed area, and K 
is a constant equal to 0.3. The potential energy is minimized by solving for the 
coordinates of the vertex, Xj and y : . Only a small number of smoothing sweeps 
are required. Extremely irregular geometries such as the trailing edges of viscous 
O-meshes can be smoothed with this method without underrelaxation. 


3.2 Grid Adaption 

A final issue to be addressed in the grid generation process is grid adaption. As 
mentioned earlier, the time required to reach a steady state on a grid of specified 
fineness can be sped up if major flow features are first resolved on a coarse grid 
and grid points are then added in the areas of high gradients only. In particular, 



CHAPTER 3. GRID GENERATION 


■lllllllll 


rngmmrnSm 



Riniffillllll 

Hna 



■OH 

■MIH 



Figure 11: Final Grid - Quadrilaterals and Triangles 
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angle of attack (degrees) 

0 

5 

10 

15 






nodes 

11034 

14880 

13972 

19595 

cells (% triangles) 

11599 (13) 

15935 (8) 

14601 (12) 

21220 (18) 

edges 

22633 

30275 

28753 

40815 

points on the airfoil 

197 

274 

281 

368 


Table 1: Grid Dimensions: h/c = 2.25, Re = 1,200,000 


boundary layers and wakes can be continually refined until a specified spacing, e.g. 
a minimum value of y + for turbulent flows, is obtained. This procedure also removes 
the problem of having to guess a priori the location of the wakes. The adaption 
program starts with a mesh and a flow solution on that mesh. The average difference 
of the refinement flow variable between endpoints of all edges is computed. If the 
undivided difference of the flow variable across an edge is greater than a specified 
percent of this average, a node is added at the midpoint of the edge [12]. The new 
points are joined in so as to form quadrilaterals whenever possible and retain any 
structure of the original mesh. Mach number and entropy are typically used as the 
refinement variables. They allow for resolution of all relevant flow features: boundary 
layers, wakes, shocks, separation, and leading edge peaks. Figure 12 shows a sequence 
of meshes generated by this refinement process. The sequence of six meshes, not all 
of which are shown, reduced the spacing at the wall by a factor of 30. The flow 
solution is not reconverged on all intermediate meshes. Grid adaption would play an 
important role in the implementation of a multigrid algorithm. 


3.3 Grids 

The final grids for the four angles of attack investigated here are shown for h/c — 2.25 
and Re = 1,200,000 in Figure 13. The dimensions of the grids are shown in Table 1. 
All the low Reynolds number grids have typical spacings at the airfoil of 6.E-05 chords 
(y + < 10) obtained by six levels of adaption. For the higher Reynolds number cases an 
extra level of adaption was required to achieve the desired wall spacing. The number 
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Figure 12: Grid Adaption 
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Figure 13: Grids for h/c — 2.25, Re = 1,200,000 
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of triangles is indicative of the amount of adaption as they are used to transition 
between levels of nested quadrilaterals. The smaller and larger height to chord ratios 
considered in this study subtract or add layers of quadrilateral cells to the wall. The 
meshes with height to chord ratio 4.5 add 1830 points (15 layers of quadrilateral 
cells) and h/c = 1.5 meshes subtract 488 points (5 layers of quadrilateral cells) from 
the h/c = 2.25 meshes. The background Cartesian wind tunnel grid has a constant 
spacing in the y direction of 0.075 chords and geometric stretching in the x direction. 
An example of the unconstrained grid is shown in Figure 14. It adds a C-mesh outer 
boundary and 1061 points to the h/c = 4.5 grid. Because the grids only change near 
the wind tunnel walls, grid effects can be ruled out as having any influence on the 
difference between solutions of varying h/c at a constant angle of attack. There are 
61 points on each wall which extend from 5 chords upstream of the airfoil leading 
edge to 15 chords downstream. 
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Figure 14: Grid for Unconstrained Flow 




Chapter 4 

Boundary Conditions 

Boundary conditions are implemented in a manner consistent with the solution method 
in the interior of the domain. This usually involves a weak variational formulation of 
the boundary conditions. A control volume is constructed at the boundary points and 
fluxes through all segments are computed. For boundary control volumes this now 
includes a boundary flux (Figure 15). In this manner boundary conditions are only 
weakly enforced unless they are explicitly specified, for example, no slip: u = v = 0. 
Several boundary conditions are required for modeling wind tunnel flows. 


4.1 Solid Walls 

On the airfoil a no slip condition is applied. In addition the wall is specified as 
adiabatic and dT/dn is set to zero in the boundary flux. The Navier-Stokes equations 
are integrated up to the wall using these specified variables. In this way continuity and 
normal momentum are automatically solved but with reduced accuracy at the wall 
due to the one-sided control volume. No extrapolations or approximate relations are 
used as in finite difference methods to enforce normal momentum and the adiabatic 

wall. 

The wind tunnel walls are modeled as solid with a slip boundary condition. Al- 
though it would be a simple matter to include the wall boundary layer by increasing 
the grid point density near the walls so as to resolve them, in this study the effects of 


32 


CHAPTER 4. BOUNDARY CONDITIONS 


33 


— boundary control volume 



Figure 15: Weak Enforcement of Boundary Conditions 
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the wall boundary layers were not considered. Possible implications of a slip bound- 
ary condition are blowing and suction at the wall to remove the boundary layer or 
canting of the walls based on momentum thickness to negate boundary layer thick- 
ness effects. The wind tunnel walls are also adiabatic. The current implementation 
of a slip wall sets the normal velocity to zero at the end of each time step and sets 
the temperature gradient in the boundary flux. The tangential velocity component 
is extrapolated from the interior. 


4.2 Inflow and Outflow 


The information required for inflow or outflow is determined by the method of char- 
acteristics. For subsonic inflow there are three incoming characteristics which neces- 
sitate the specification of three quantities. The specified variables are the velocity 
tangential to the boundary, entropy, and the incoming Riemann invariant. The tan- 
gential velocity depends on the orientation of the boundary and is computed using 
u = Uqo, v = 0. For the wind tunnel inflow plane, this is a parallel inflow assumption, 
^tangential = v = 0, and is valid if the inflow plane is located far enough upstream. 
Entropy is set to its freestream value. The incoming Riemann invariant, 

Ri = (17) 

7-1 


is the information carried on the one-dimensional, isentropic incoming characteris- 
tic [13], where u nao is the velocity normal to the boundary at infinity. For the parallel 
inflow assumption, this equates to freestream. The outgoing Riemann invariant, 


#2 = 


2c e *t r 

7-1 


( 18 ) 


is extrapolated from the interior. The normal velocity and speed of sound at the 
boundary are then determined from 


u n 


\(Ri + R2) 


7-1 


(i?2 - Ri) 


C 


2 


( 19 ) 
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The Cartesian velocities can be determined from the components of velocity nor- 
mal and tangential to the boundary. Pressure and density are computed from the 
definition of the speed of sound c 2 = 7 p/p and the isentropic relation s = p/p 7 . 

Exit boundary conditions for internal flow problems can be difficult to formulate 
unless outside information is available. For subsonic outflow one quantity must be 
specified due to the one incoming characteristic. Extrapolation of all quantities is 
unphysical and can prevent convergence or give poor results which depend on the 
initial conditions [14, 15]. Ideally, experimental wind tunnel pressures are desired. 
Since the calculations are not being directly compared to experimental results, there 
is no obvious quantity to specify at the outflow boundary. Therefore, the boundary 
was moved far downstream (15 chords) to minimize its influence on the rest of the 
flow field, and pressure was set to its freestream value. To damp out oscillations, a 
non-reflecting boundary condition based on the incoming characteristic is used [15]: 

dp du , . 

- pc-fo + <*(p - Poo) = 0 (20) 

or numerically: 

P" +1 = P] + P n c"A*K +1 - «?) - «At(p? - Poo ) (21) 

This helps but does not eliminate reflection of pressure disturbances back into the 
wind tunnel where they bounce between the exit and entrance plane and the walls. At 
steady state ( d/dt = 0) the outflow pressure reaches The remaining variables are 
extrapolated. The inflow and outflow boundary conditions used for the wind tunnel 
case can also be used in the unconstrained case. For the C-mesh used on the outer 
boundary, the orientation of the inflow boundary now determines the tangential and 
normal directions which are not always aligned with the Cartesian directions. 



Chapter 5 
Discussion 

Four angles of attack (a = 0, 5, 10, and 15 degrees), three total wind tunnel 
height to airfoil chord ratios ( h/c = 1.5, 2.25, and 4.5), and two Reynolds num- 
bers (Re = 1,200,000 and 6,000,000) were computed. The intermediate height to 
chord ratio of 2.25 was not run for the higher Reynolds number. The unconstrained 
case was simulated using an outer boundary diameter of 24 chords and appropriate 
boundary conditions. All cases use a freestream Mach number of 0.25. In this chapter, 
the results from the CFD runs will be presented and discussed. 


5.1 Results 

At zero degrees angle of attack, the airfoil is generating a small amount of lift. Fig- 
ure 16 shows negative pressure coefficients on the airfoil for h/c = 1.5, 2.25, and the 
unconstrained case at the lower Reynolds number while Figure 17 compares distri- 
butions at the two Reynolds numbers. The pressure distributions vary only slightly. 
The stagnation point at zero degrees angle of attack is right at the leading edge 
and bumps and flat spots in the geometry are indicated by the bumps and spikes in 
the solution on both the pressure and suction sides of the airfoil. The pressure side 
leading edge spike is seen to go away as the stagnation point moves back at higher 
angles of attack. Plots of skin friction coefficient, C/, are shown in Figures 18 and 19 
for the low and high Reynolds numbers respectively. Trends are similar, but the low 
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Figure 20: Wall Pressures: a = 0 degrees 


Reynolds number case shows a small amount of trailing edge separation and increased 
Cf. Portions of the upper and lower surfaces would be laminar if a transition model 
were implemented, i.e. if the flow were not tripped at the leading edge. Figure 20 is 
a plot of wall pressures for the confined flow cases. The effects here can be considered 
as the addition of a potential source and vortex. The source represents blockage. The 
vortex adds the small amount of lift. Effects from a source would be symmetric about 
the centerline on the top and bottom walls, while the vortex would show an antisym- 
metric distribution. Both are seen to influence the wall distribution. Because of the 
Riemann invariant boundary conditions at the inflow, the incoming Mach number 
and pressure will not necessarily be freestream values, but the exit pressure decays 
to the prescribed freestream value. 

Similar plots of airfoil pressure distributions, skin friction coefficients, and wall 
pressure distributions for the five degrees angle of attack case are shown in Figures 21 
through 25. As expected the lift increases as the walls are moved in (Figure 21), but 
Reynolds number is seen to have very little effect on the airfoil pressure distributions 
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(Figure 22). The skin friction plots for the two Reynolds numbers, Figures 23 
and 24, indicate no effect of wall proximity on the small recirculation region present 
at the trailing edge except for the unconstrained, low Reynolds number case where 
this region is slightly reduced. The recirculation region begins at the separation point 
where C f equals zero and extends to the trailing edge. The region is indicated by the 
negative, almost zero skin friction coefficient and is smaller at the higher Reynolds 
number. The major effect seen in the wall pressures is that of a vortex (Figure 25). 
Reynolds number has no influence at the slip walls. Typical flowfields are shown in 
the Mach number contour plots of Figures 26 and 27 at five degrees angle of attack 
for the the height to chord ratio of 1.5 and the unconstrained runs. Basic qualitative 
features of the wind tunnel flows are illustrated, for example, reduced wake curvature 
and the confining effects of the walls. 

Figure 28 shows the effect of wall proximity for the ten degrees angle of attack case 
and Re = 1,200,000. Effects of higher Reynolds number are shown in Figure 29 where 
some differences are apparent for the constrained cases, especially the lowest height 
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to chord ratio. The grid adaption allows the leading edge peak to be accurately 
captured. Figures 30 and 31 show skin friction plot for Reynolds numbers of 1.2 
and 6.0 million, respectively. They indicate a reduced recirculation region at the 
trailing edge as h/c is increased. This is due to the decreased levels of lift which delay 
separation. Wall pressures are shown in Figure 32 where the Reynolds number effects 
at the airfoil are reflected. 

Results for the fifteen degrees angle of attack case are presented in Figures 33 
through 38. Figures 33 and 34 are airfoil C p distributions. The critical pressure 
coefficient, C*, at M = 1.0 is —10.2, so that for this configuration the flow becomes 
supersonic over the leading edge. Figure 35, a Mach number contour plot of the 
leading edge of the airfoil for the high Reynolds number, small height to chord ratio 
case, shows a small shock which appears at this angle of attack. This is probably 
caused in part by the irregular geometry in this region. The boundary layer thickens 
considerably behind the shock. For h/c — 1.5 and Re = 6,000,000 the peak Mach 
number is 1.70. For Re = 1,200,000 the peak Mach number is 1.62, while for the 
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Figure 24: Airfoil Skin Friction: a = 5 degrees, Re = 6,000,000 
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Figure 25: Wall Pressures: a = 5 degrees 


unconstrained case at this Reynolds number it is 1.45. The skin friction plot for the 
lower Reynolds number (Figure 36) indicates leading edge separation, reattachment, 
and separation again over much of the rear of the airfoil with the extent of the 
separation dependent on height to chord ratio. For the higher Reynolds number 
the leading edge separation is reduced (Figure 37). Wall pressures are shown in 
Figure 38 for both Reynolds numbers. Despite the differences in the behavior of 
the skin friction, Reynolds number has relatively little effect on the airfoil and wall 
pressure distributions. The wake curvature is more apparent in the Mach number 
contour plots of Figures 39 and 40 with the airfoil at fifteen degrees angle of attack. 
Here the larger separation region and thicker boundary layer and wake of the lowest 
height to chord ratio configuration are compared with the unconstrained solution. 




CHAPTER 5. DISCUSSION 


44 



Figure 27: Mach Number Contours: a = 5 degrees, unconstrained, Re = 1,200,000 
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5.2 Lift 

Table 2 summarizes the lift and drag coefficients as a function of wind tunnel height, 
Reynolds number, and angle of attack. Trends in the lift coefficient show that in all 
cases it increases with decreasing wind tunnel wall height to chord ratio. This is due 
to the wall induced upwash [16] and the vortex reflection and acceleration effects the 
walls have on the flow. Increasing the Reynolds number tended to increase the lift 
coefficient slightly, by at most 6%, with the largest percentage changes occurring at 
small h/c. The small variations of pressure coefficient with Reynolds number are seen 
in Figures 17, 22, 29, and 34. It seems that increasing the Reynolds number has, to 
a much lesser extent, the same effect on lift as decreasing the height to chord ratio. 

Figure 41 shows lift coefficient versus angle of attack for the four height to chord 
ratios. Up to ten degrees where the flow is almost fully attached, the data indicates 
the effect of decreasing h/c or increasing Reynolds number is to linearly displace 
the lift curve upward and increase its slope. Even in the region near stall past ten 
degrees angle of attack where the lift curves level off and viscous phenomenon play 
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Figure 35: Blow Up of Leading Edge Mach Number Contours: a = 15 degrees, 
h/c = 1.5, Re = 6,000,000 
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Table 2: Lift and Drag Coefficients 
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Figure 41: Lift Curves versus h/c 
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Table 3: Pressure and Viscous Drag Coefficients 


an important role in determining the flow field, there is no indication of any signif- 
icant nonlinear behavior of lift coefficient with changes in h/c. Wind tunnel wall 
interference correction methods which are based on the idea of global corrections to 
freestream Mach number and angle of attack seem to be supported since lift curve 
slope is a function of Mach number in two dimensions and lift level is a function of 
angle of attack. But while the idea of AAf and A at may be correct, the effects which 
go into calculating these corrections need to be considered. 

At the walls, increasing the Reynolds number shows no discernible trends (Figures 
20, 25, 32, and 38) and in many cases has no effect at all on wall pressures, especially 
at the highest height to chord ratio. Figure 41 shows there to be some Reynolds 
number effects at the airfoil for h/c = 4.5 even when wall pressures show none. 

5.3 Drag 

The drag can be broken down into that due to viscous forces, Cd, itc , and pressure 
forces, Cd ^.,.1 and these values are indicated in Table 3. Viscous drag is due to the 
viscous shearing forces and acts tangential to the solid surface. Pressure or form drag 
acts normal to a surface and is due to the displacement of the streamlines in the 
presence of separation. It is seen that the pressure drag makes a large contribution to 
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the total drag coefficient when the airfoil is at high angles of attack and is no longer 
a streamlined body. 

Two major effects account for the variations in viscous drag behavior. The first 
effect causes an increase in with decreasing height to chord ratio. The increase 

is due to the higher speed flow over the airfoil which increases the boundary layer edge 
velocity and at the same time the wall shear and, therefore, the skin friction. This 
is the dominant effect seen in all of the lower angle of attack, attached, constrained 
cases, where the viscous drag increases as the height to chord ratio is reduced (Figures 
18, 19, 23, 24, 30, and 31). 

The second, competing effect contributes to a decreased value of viscous drag with 
decreasing height to chord ratio. Since the velocities over the airfoil are increased as 
the height to chord ratio is decreased, the higher levels of lift promote separation 
which enlarges the recirculation region at the trailing edge. Increased amounts of 
reversed flow in the separated region decrease viscous drag by adding what can be 
considered as thrust. At fifteen degrees angle of attack the amount of separation 
is significant and the point of separation changes considerably with hjc (Figures 36 
and 37). For the high angle of attack at both Reynolds numbers, this is the major 
influence on viscous drag, causing it to decrease with decreasing height to chord ratio. 
This is also the main effect for the attached, unconstrained, low Reynolds number 
runs, where the viscous drag trend reverses from that of the constrained cases and 
Cd* itc increases because the amount of separation is reduced. For the high Reynolds 
number cases the trend does not reverse and the viscous drag further decreases in 
going from the constrained to unconstrained case. 

In general at both Reynolds numbers, pressure drag follows the same trends as 
the viscous drag except at fifteen degrees angle of attack where it shows the oppo- 
site trend. Blockage effects of the boundary layer and larger recirculation regions 
substantially increase pressure drag so that form drag makes the overall drag show a 
marked increase as hjc is decreased and the separation point is moved upstream on 
the airfoil. In all cases the trailing edge separation is delayed at the higher Reynolds 
number. 
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The influences on drag are then the increased velocity over the airfoil which in- 
creases lift and viscous drag but also encourages separation which decreases lift and 
viscous drag but substantially increases pressure drag. Total drag coefficients gener- 
ally increase with decreasing h/c, but are lowest for the highest height to chord ratio 
rather than the unconstrained case. 


5.4 Code Characteristics 

It should be noted that convergence to steady state was very slow for all cases. 
Because the code is explicit and high aspect ratio boundary layer cells make the 
equations stiff, thousands of iterations are needed to obtain a converged solution. 
Typically, 35,000 iterations were required on a sequence of grids for approximately 
four orders of magnitude reduction in the residual and a lift coefficient converged 
to within 1%. Total CPU time per run is approximately seven hours, depending on 
the size of the grid, using a single processor of a Cray YMP. CFL numbers based 
on monotonicity principles as high as 2.3 were used. The entire code is vectorized 
for increased speed on the Cray by “coloring” edges and faces so that loops with 
data dependencies are broken down into several nonrecursive loops. Solutions started 
from other solutions, e.g. when the height to chord ratio was reduced or the Reynolds 
number was increased, required approximately 10,000—15,000 iterations on the finest 
grid for reconvergence. 

The convergence is especially slow for two reasons. First, the formulation is com- 
pressible and all of the solutions presented here are for small Mach number. Current 
algorithms developed for compressible flow based on wave propagation do not effi- 
ciently model the elliptic nature of low Mach number flow. The confined wind tunnel 
also contributed to convergence problems. Pressure waves reflect inside the tunnel, 
and even when they are damped with appropriate boundary conditions, they slow 
convergence and cause the solution to oscillate. Low frequency waves are the hardest 
to d am p out and their presence is not indicated by the residual which only measures 
local, high frequency error. Low frequency pressure oscillations can be seen in Fig- 
ure 32 on the walls behind the airfoil for h/c = 1.5. Although lift and viscous drag 
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were converged to within 1% for all cases, small changes in surface pressures produce 
large changes in the pressure drag. For the smaller height to chord ratios at higher 
angles of attack, the solutions often never did converge enough to give accurate values 
of pressure drag. This difficulty could be attributed to the surface pressure integra- 
tions, but the surface grids are fairly refined, especially at the leading edge due to the 
grid adaption. Similar convergence problems for nearly incompressible, confined flow 
have been reported on structured grids using implicit, compressible methods [17]. The 
unconstrained flow cases showed more monotonic convergence as did the supercritical 
cases. 

The main reason for using the present code was its ability model complex geome- 
tries in anticipation of computing the flow about multielement airfoils. In its current 
form, though, it can only be considered a research tool. 



Chapter 6 
Conclusions 


A method for computing wind tunnel flows has been developed using an unstructured 
grid Navier-Stokes computational fluid dynamics code. The code uses the Baldwin- 
Lomax algebraic turbulence model. For the unstructured grid, boundary layer profiles 
normal to the body are constructed using boundary geometry information. The grid 
generation process has been illustrated on a complex multielement airfoil geometry to 
show the benefits of an unstructured grid while still retaining some of the advantages 
offered by structured grids, namely the stretched quadrilateral boundary layer and 
wake cells. Special procedures are employed to handle these high aspect ratio cells 
during the triangulation process. Grid adaption was used to resolve various flow 
features and to speed up convergence to steady state. Using grid adaption on an initial 
coarse grid helps ease some of the difficulties of working with the highly stretched 
cells. The exit boundary condition is influential in wind tunnel flow computations, 
and without the advantage of experimental wind tunnel test data, freestream pressure 
was imposed at the outflow boundary. 

The ability to model wind tunnel flows directly is useful for CFD validation, 
particularly turbulence models; wind tunnel wall correction methods development; 
and fundamental studies of flow physics. Results were presented of an investigation 
of the effect of wind tunnel wall proximity and Reynolds number on a single element 
airfoil model at low speed and four angles of attack up to stall. Decreasing the wind 
tunnel height to chord ratio was seen to increase both the lift and drag coefficients at 
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all angles of attack, although the drag coefficient also increased for the unconstrained 
cases from that of the h/c = 4.5 cases. Closer wind tunnel wall proximity and lower 
Reynolds number promoted earlier separation. Separation decreases viscous drag, 
but the larger extent of the boundary layer and recirculation region significantly 
increase pressure drag which then becomes the major contributor to total drag. For 
the subsonic cases considered, the results support the idea of Mach number and angle 
of attack corrections to wind tunnel data to obtain an equivalent unconstrained flow. 
For the attached flow cases, the lift curves remained linear but their slopes and levels 
of lift increased as the wind tunnel walls were moved in or as Reynolds number was 
increased. 

The highest height to chord ratio of 4.5 does not simulate the unconstrained flow 
at any of the angles of attack investigated here as indicated by the airfoil pressure 
distributions or global aerodynamic parameters. Wind tunnel wail interference cor- 
rection methods would be required for ail cases. The effects which go into computing 
the corrections should include viscous phenomena for optimal accuracy, especially at 
higher angles of attack. Wall pressure plots did not always reflect Reynolds number 
effects at the airfoil. It seems doubtful that wind tunnel wall interference correction 
methods which use only wall pressure measurements could differentiate Reynolds 
number effects at the airfoil when no viscous modeling is included. Correction meth- 
ods should therefore use airfoil pressures as well as wall measurements or include a 
viscous formulation in the numerical method. 

Future work should focus on enhancements to speed up the code. Computation 
times for unstructured grid, explicit codes with multigrid have been demonstrated 
to approach the speed of equivalent structured, implicit codes [7]. The benefits of 
unstructured grids, including the ability to model complex flows and grid adaption to 
relevant flow features, make some additional time penalty worthwhile. Modifications 
to the Baldwin-Lomax turbulence model to include the ability to handle confluent 
boundary layers or a switch to a one or two equation model is also recommended. 
Changes to the algebraic model require being able to differentiate between boundary 
layers and overlying wakes so that separate length scales can be computed. Devel- 
opment of a porous wall boundary condition will allow investigation of porous wall 
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wind tunnels, although they are rarely used for the low Mach number flows considered 
here. Resolution of the wall boundary layer will distinguish the airfoil viscous effects 
from those due to the walls. 
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